Active contours and active shape models (ASM) have been widely employed in image segmentation. ASMs are, however, limited in their inability to (a) resolve boundaries of intersecting objects and (b) to handle occlusion. Multiple overlapping objects are typically segmented out as a single object. In addition, ASMs are limited by point correspondence issues, since object landmarks need to be identified across multiple objects for initial object alignment. ASMs also are constrained in that they can usually only segment a single object in an image and hence are unable to resolve multiple overlapping objects.
Statistical Shape Models (SSM) (Cootes, T. F. et al., “Active Shape Models—Their Training and Application,” Computer Vision and Image Understanding, January 1995, 61(1): 38-59) and Active Contours (AC) have been used heavily to segment objects of interest in an array of applications involving biomedical imaging. (Veltri, R. W. et al., “Nuclear Roundness Variance Predicts Prostate Cancer Progression, Metastasis, and Death: A Prospective Evaluation With up to 25 Years of Follow-Up After Radical Prostatectomy,” The Prostate, 2010, 70, 1333-1339; Fatakdawala, H. et al., “Expectation Maximization driven Geodesic Active Contour with Overlap Resolution (EMaGACOR): Application to Lymphocyte Segmentation on Breast Cancer Histopathology,” Biomedical Engineering, IEEE Transactions on, 2010, 57(7): 1676-1689).
In recent years, shape based active contour models have emerged as a way to overcome overlap resolution. However, most of these shape-based methods are limited to finding and resolving one object overlap per scene and require user intervention for model initialization. A fundamental problem when computing SSMs is the determination of correspondences between training instances of the object of interest, especially when the surfaces are represented by point clouds. (Cootes, T. F. et al., “Active Shape Models—Their Training and Application,” Computer Vision and Image Understanding, January 1995, 61(1): 38-59). Often, homologies between the points are assumed, which might lead to an imprecise model of the shape.
A number of deformable segmentation schemes have been developed to date. They can be divided into two categories: boundary-based (first generation) and region-based (second generation) models.
Boundary-based approaches such as geodesic/geometric active contours (ACs) are popular for their reliable performance and fast computation. (Kass, M. et al., “Snakes: active contour models,” International Journal of Computer Vision, 1987, pp. 321-331; Caselles, V. et al., “Geodesic active contours,” Int. J. Comput. Vision, 1997, 22(1): 61-79; Kichenassamy, S. wt al., “Conformal curvature flows: From phase transitions to active vision,” In Archive for Rational Mechanics and Analysis, 1996, 134(3): 275-301, 1996). However, as only the edge information is utilized, their performance is limited by the strength of the image gradient. Such boundary-based models are typically unable to handle object occlusion or scene clutter.
AC models based on edge/boundary detection are typically referred to as the first generation of active contours. (Kass, M. et al., “Snakes: active contour models,” International Journal of Computer Vision, 1987, pp. 321-331; Caselles, V. et al., “Geodesic active contours,” Int. J. Comput. Vision, 1997, 22(1): 61-79).
Kass et al. described a boundary-based geodesic active contour model. According to the Kass et al. model, the segmentation of any object in a given image, which is well discernible and whose edges can be described by a closed curve, is equivalent to the location of sharp image intensity variations by iteratively deforming a curve  towards the edges of the object. Such a model is entirely dependent on the chosen parametrization of the initial curve (p)=(x(p), y(p))εΩ, pε[0,1]
Caselles et al. described a method to solve the parametrization dependency problem with the Geodesic Active Contour (GAC) model. Caselles et al. introduced a new intrinsic energy functional which is independent of the initial curve parametrization, and minimization of which results in a geodesic in Riemannian space. The associated evolution equation can be written as ∂tC=F·, where F is a speed term derived from an energy functional and  is the exterior unit normal vector to the curve . ∂t is the partial derivative at an instance of time t. C represents a two-dimensional Cartesian grid of pixels c=(x, y).
To overcome limitations associated with topological changes of object boundaries, such as splitting and merging, the evolution curve is handled by a level set method described by Osher et al. (Osher, S. et al., “Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations,” Journal of Computational Physics, 1988, 79(1): 12-49). The level set method includes representing the curve as the level zero of a higher dimensional function φ, which is often a signed distance function. Thus, each active contour (AC) or active shape within a given training set of sample shapes is embedded as a zero level set of a higher dimensional surface. A Signed Distance Function (SDF) is a functional that depends on the AC providing the boundaries. The level set function φ evolves in time as: ∂tφ=F·|Δφ| and the interface at the level zero corresponds to the model's evolving contour.
However, the use of edge-based active contours in segmentation methods is limited because of high sensitivity to noise, requiring careful initialization. To overcome such limitations, segmentation models based on region statistics such as mean, variance, and probability density functions (PDFs) have been developed.
Chan (2005) proposed a 2-phase segmentation method based on a mean descriptor where the AC evolves in such a way that the difference between the gray level intensity average between the foreground and the background was maximized. (Chan, T. “Level set based shape prior segmentation,” IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2005, 2: 1164-1170). Region-based approaches rely on statistics obtained from the entire region. Color or textural attributes being often employed to detect the object. (Chan, T. F. et al., “Active contours without edges,” IEEE Trans. on Image Processing, 2001, 10(2): 266-277). However, these models typically require far more computations and, like the boundary-based approaches, are limited in their ability to deal with object occlusions and scene clutter.
Third generation AC models involve combining a shape prior with geometric/geodesic active contours that simultaneously achieves registration and segmentation. A number of methods incorporating various shape priors into active contour formulation to resolve overlaps have been described.
Cootes et al. (1995) described a principal component analysis (PCA) method to capture the main shape variations of parametric active contours (active shapes) within a training set and thus to represent the prior shape information. (Cootes, T. F. et al., “Active Shape Models—Their Training and Application,” Computer Vision and Image Understanding, January 1995, 61(1): 38-59). Consequently, the Cootes et al. model is not parametrization free.
Leventon, et al. described a method of introducing prior shape information into AC, intrinsically represented by level set functions. Such a method applies PCA on the signed distance functions (SDF) of the parametric Geodesic Active Contours (GACs). (Leventon et al., “Statistical shape influence in geodesic active contours,” in CVJR, 2000, pp. 316-323). Such a feature allows construction of an intrinsic and parametrization free shape model. SDFs have the additional advantage that they are more robust to slight misalignments of the training sequence compared to parametric curves. However, the shape functions resulting from a PCA are not exactly SDFs, but are nonetheless used in practice since they are close to real SDFs. Rousson et al (2002) described a method for introducing shape priors into level set representations targeting 2D closed structures. (Rousson, M. et al., “Shape priors for level set representation,” European Conference on Computer Vision, 2002, pp. 78-92).
Chan described a method of combining a statistical shape prior with geometric/geodesic active contours via a variational model that allows for simultaneous registration and segmentation. (Chan, T., “Level set based shape prior segmentation”, IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2005, 2: 1164-1170). Bresson, et al. (Bresson, X. et al., “A priori information in image segmentation: energy functional based on shape statistical model and image information,” IEEE Int. Conf. on Image Processing, September 2003, 3: 425-428) described a method integrating a geometric shape prior described in Leventon et al. (Leventon et al., “Statistical shape influence in geodesic active contours,” in CVJR, 2000, pp. 316-323) into the segmentation framework based on AC as well as on a region driven energy term described in Chan (Chan, T. “Level set based shape prior segmentation,” IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2005, 2: 1164-1170). Paragios et al. (Paragios, N. et al., “Unifying boundary and region-based information for geodesic active tracking,” IEEE Conference on Computer Vision and Pattern Recognition, 1999, 2: 300-305) and Rousson et al (Rousson, M. et al., “Shape priors for level set representation,” European Conference on Computer Vision, 2002, pp. 78-92) describe a method of computing signed distance functions of the training images and capturing the statistics of the signed distance training set via the PCA method. Such a method assumes that the probability distribution function (PDF) of the training set is Gaussian.
A limitation of third generation AC models is that they introduce shape priors into a level set framework in such a way that usually only one pair of overlapping objects can be accurately resolved into independent shapes within an image scene. (Fang, W. et al., “Statistical Shape Influence in Geodesic Active Contours,” IEEE Conference on Computer Vision and Pattern, 2007, 40(8): 2163-2172; Bresson, X. et al., “A priori information in image segmentation: energy functional based on shape statistical model and image information,” IEEE Int. Conf. on Image Processing, September 2003, 3: 25-428) Further, all third generation AC methods are sensitive to model initialization and typically require varying degrees of user intervention.
The flexibility of the level set method has traditionally resulted in long computation times (up to several hours for large images with several thousand structures to segment) and therefore of limited practical utility in dealing with very large histopathological images. As a result, processing times in a clinical setting may be impractical. Thus there is a clear need to accelerate these segmentation schemes by using a parallel algorithm executed on a Graphical Processing Unit (GPU). Lefohn et al. (Lefohn A. E. et al., “GPU based, three-dimensional level set solver with curvature flow,” Technical report uucs-02-017, School of Computing, University of Utah, (2002)) and Tejada et al. (Tejada E. et al., “Large steps in GPU-based deformable bodies simulation,” Simulat Model Pract Theory, 2005, 13(8):703-715) describe a method using graphical processing unit (GPU) architectures employing higher-level segmentation approaches, such as deformable models, by considering implicit deformable models such as image lattices (e.g. a 2D curve is implicitly represented as the iso-value of a field encoded as a 2D image). Level-sets approaches became particularly popular in the GPU-segmentation community as significant speed-ups and interactive rendering were made available. Geodesic active contours, which are a combination of traditional active contours (snakes) and level-sets evolution, were efficiently implemented in GPU by using the total variation formulation, mostly to quickly segment structures by foreground and background separation in 2D images. (Schmid, J. et al., “A GPU framework for parallel segmentation of volumetric images using discrete deformable model,” The Visual Computer, 2011, 27(2): 85-95). Ruiz et al. described a GPU accelerated segmentation schemes with applications to biomedical histopathology. (Ruiz, A. et al., “Pathological image segmentation for neuroblastoma using the GPU,” Biomedical Imaging: From Nano to Macro, IEEE International Symposium on, pp. 296-299, (2008).)
In the rapidly emerging field of digital pathology, the ability to segment multiple objects, especially the ability to deal with object overlap and occlusion, is highly critical in the context of a number of different diagnostic and prognostic applications. (Gurcan, M. et al., “Histopathological Image Analysis: A Review,” IEEE Reviews in Biomedical Engineering, 2009, 2: 147-171; Madabhushi, A, “Digital Pathology Image Analysis: Opportunities and Challenges,” (Editorial), Imaging in Medicine, 2009, 1(1): 7-10).
Currently, the diagnosis of diseases, such as prostate cancer and breast cancer, is done manually by visual analysis of tissue samples, typically obtained from a patient via biopsy. (Fatakdawala, H. et al., “Expectation Maximization driven Geodesic Active Contour with Overlap Resolution (EMaGACOR): Application to Lymphocyte Segmentation on Breast Cancer Histopathology,” Biomedical Engineering, IEEE Transactions on, 2010, 57(7): 1676-1689; Basavanhally, A. et al., “Computerized Image-Based Detection and Grading of Lymphocytic Infiltration in HER2+ Breast Cancer Histopathology,” Biomedical Engineering, IEEE Transactions on, 2010, 57 (3): 642-653). The architectural arrangement of nuclear and glandular structures on histopathology is highly relevant in the context of disease grading. (Fatakdawala, H. et al., “Expectation Maximization driven Geodesic Active Contour with Overlap Resolution (EMaGACOR): Application to Lymphocyte Segmentation on Breast Cancer Histopathology,” Biomedical Engineering, IEEE Transactions on, 2010, 57(7): 1676-1689; Doyle, S. et al., “Automated Grading of Prostate Cancer using Architectural and Textural Image Features,” International Symposium on Biomedical Imaging (ISBI), 2007, pp. 1284-87; Basavanhally, A. et al., “Computerized Image-Based Detection and Grading of Lymphocytic Infiltration in HER2+ Breast Cancer Histopathology,” Biomedical Engineering, IEEE Transactions on, 2010, 57 (3): 642-653).
Cancer grade in the context of breast and prostate cancer is a key feature used to predict patient prognosis and in prescribing a treatment. (Gurcan, M. et al., “Histopathological Image Analysis: A Review,” IEEE Reviews in Biomedical Engineering, 2009, 2: 147-171). While grading systems exist in the context of both breast cancer and prostate cancer, manual grading is time consuming and prone to human errors due to observer variability and can lead to variable prognosis and suboptimal treatment. Automated segmentation and quantification of nuclear and glandular structures is critical for classification and grading of cancer. (Fatakdawala, H. et al., “Expectation Maximization driven Geodesic Active Contour with Overlap Resolution (EMaGACOR): Application to Lymphocyte Segmentation on Breast Cancer Histopathology,” Biomedical Engineering, IEEE Transactions on, 2010, 57(7): 1676-1689; Doyle, S. et al., “Automated Grading of Prostate Cancer using Architectural and Textural Image Features,” International Symposium on Biomedical Imaging (ISBI), 2007, pp. 1284-87).
In the context of prostate cancer (CaP), pathologists grade histopathological specimens by visually characterizing gland morphology and architecture in regions they suspect are malignant. The Gleason grading system is used to describe CaP aggressiveness; lower Gleason grade structures, such as glands, are medium-sized with round shapes, while in higher Gleason grade patterns, such as glands, tend to be small and have irregular shapes. Gleason grading ranges from very well differentiated (grade 1) to very poorly differentiated (grade 5). Doyle, et al. (2007) showed that spatial graphs (eg. Voronoi, Delaunay, minimum spanning tree) built using nuclei as vertices in digitized histopathology images yielded a set of quantitative feature that allowed for improved separation between intermediate gleason patterns. (Doyle, S. et al., “Automated Grading of Prostate Cancer using Architectural and Textural. Image Features,” International Symposium on Biomedical Imaging (ISBI), 2007, pp. 1284-87). Veltri, et al. (2010) showed that nuclear shape and morphology was reflective of disease aggressiveness and patient outcome. (Veltri, R. W. et al., “Nuclear Roundness Variance Predicts Prostate Cancer Progression, Metastasis, and Death: A Prospective Evaluation With up to 25 Years of Follow-Up After Radical Prostatectomy,” The Prostate, 2010, 70, 1333-1339). Both of these methods require accurate segmentation of prostate nuclei as an initial step, but had previously employed manual or semi-automated approaches.
Lymphocyte Infiltration (LI) has been identified as an important prognostic marker of outcome in Her2+ breast cancer and in other diseases as well. (Basavanhally, A. et al., “Computerized Image-Based Detection and Grading of Lymphocytic Infiltration in HER2+ Breast Cancer Histopathology,” Biomedical Engineering, IEEE Transactions on, 2010, 57 (3): 642-653; Zhang, L. et al., “Intratumoral T Cells, Recurrence, and Survival in Epithelial Ovarian Cancer,” New England. Journal of Medicine, 2003, 348(3): 203-213). Lymphocyte segmentation in hematoxylin and eosin (H&E) stained BC histopathology images is complicated by the similarity in appearance between lymphocyte nuclei and other structures (e.g., cancer nuclei) in the image. This may lead to incorrect determination of the extent of LI. Automated detection and quantification of these structures on histopathology imagery could potentially result in the development of a digital prognostic tool for predicting patient outcome (and hence deciding treatment strategies). However, an automated lymphocyte or nuclear detection algorithm on H&E images has to be able to deal with (1) variability in digital slide appearance due to inconsistencies in histological staining, (2) poor image quality with tissue samples due to slide digitization, and (3) tissue fixation. Moreover, LI may be characterized by a high density of lymphocytes or nuclei, which can cause significant overlap among lymphocyte nuclei and other structures on the H&E images. Such a method often results in adjacent nuclei visually appearing as one single lymphocyte. Basavanhally, et al. quantified the extent of lymphocytic infiltration (LI) in HER2+ breast cancers using a nuclear detection and graph feature based approach. (Basavanhally, A. et al., “Computerized Image-Based Detection and Grading of Lymphocytic Infiltration in HER2+ Breast Cancer Histopathology,” Biomedical Engineering, IEEE Transactions on, 2010, 57 (3): 642-653). Fatakdwala, et al. combined an expectation maximization scheme with an explicit concavity based overlap resolution scheme to separate overlapping lymphocytic nuclei. (Fatakdawala, H. et al., “Expectation Maximization driven Geodesic Active Contour with Overlap Resolution (EMaGACOR): Application to Lymphocyte Segmentation on Breast Cancer Histopathology,” Biomedical Engineering, IEEE Transactions on, 2010, 57(7): 1676-1689).
The present invention provides a method for segmenting single non-overlapping objects or multiple overlapping/occluded objects in an image scene via a hybrid adaptive active contour model. The method synergistically combines shape priors with boundary and region-based active contours in a level set formulation with a watershed scheme for model initialization for identifying and resolving multiple object overlaps in an image scene. The method comprises learning a shape prior for a target object of interest and applying the shape prior into a hybrid active contour model to segment overlapping and non-overlapping objects simultaneously.
The hybrid active contour model is able to handle overlaps between multiple intersecting and adjacent objects and is implemented through a graphical processing unit (GPU) framework comprises representing active contours of an object in an image scene using the level set method. The level set method comprises representing the contour implicitly as the zero level of a higher dimensional embedding function, and performing contour propagation by evolving the embedding function by minimizing variational energy. The method allows ease of handling of topological changes of object boundary such as splitting or merging. By leveraging the commodity GPU, the segmentation method provides real-time speed ups, which are otherwise masked by latencies incurred due to the complexities of the algorithm.